source ReplicaSim.tcl

#setup Nodes
ReplicaSim create TERMD
TERMD setup { bin /home/yi/simulation/software/mbuilder/bin/imm1pb}
TERMD hasDimension temp
TERMD hasDimension com
TERMD setDimension temp { 298 318 328}
TERMD setDimension com { 12.0 16.0 20.0}
TERMD addRelatedVarTo com force {10 10 10 10 10}
TERMD addRelatedVarTo com win [ getWinWidth { 10 12 14 16 18}]
#startup Nodes
TERMD startNodes
TERMD do {
    load rtf ~/simulation/common/toph19_eef1.1.gouy.inp
    load para ~/simulation/common/param19_eef1.1.gouy.inp
    load psf pept.psf
}
#configure solvent,TERMDstep
IMM1_memb configure { thick 27}
MD create TERMDstep {
    numsteps 200
    timestep 2
    initTemp 298
    deltaTemp  10.0
    heatFreq 1000
    temperature 298
    COMmotion yes
    seed [randomseed]
    langevin on
    langevinFriction 5.0
    nbondedFreq 20
    rescaleFreq 10
    saveFreq  100
    printFreq 100
    restartFile run.res
    dcdFile run.dcd    
}
# create com constraint
GeoCons create slice 
TERMD copyObj { IMM1_memb TERMDstep slice}
#start simulation
TERMD do {
   setSolvent IMM1_memb {
      temperature 298
   }
   TERMDstep configure {temperature $temp seed [randomseed]}
   load coor [randompick {o1.crd o2.crd o3.crd o4.crd o5.crd o6.crd}]

   #run TERMDstep {temperature $temp}
}

TERMD do {
    move2 all -z $com
    turnon slice -z $com force $force window $win

}
# just to label that you have already finished setups and ready to do simulations.
# Used to help Replica analysis identify the setup codes above.
TERMD begin equlibrate

#TERMD do { minimize sd 300 } until { $ener(GMRS) < 1}
#   sendCmd "MMFP\nGEO plane -\n zref $com zdir 1.0 force $force droff $win sele all end\n END\n"

# Each node run simulation with constant variables but exchanging conformation
# Use as ReplicaSim run ?numsteps ?StepSetup exchange ?dimension-to-be-exchanged when { ?condition}
# Currently only one dimension can be swapping.

# equlibrate steps
TERMD run 3 TERMDstep exchange temp when {
set k 1.98e-3
puts "exchanging $node1\($temp($node1)\)-$node2\($temp($node2)\)?"

set delta [expr $energy($node1,ENER)-$energy($node2,ENER)]
set delta [expr $delta/(1/($k*$temp($node2)- 1/($k*$temp($node1))))]
set probability [expr exp(-$delta)]
puts "$energy($node1,ENER) $energy($node2,ENER) \nprobability to exchange [expr exp(-$delta)]"

if {$probability >= 1 } { expr 1 } else { expr [prob $probability]  }
}

TERMD run 5 TERMDstep exchange temp when {
set k 1.98e-3
puts "exchanging $node1\($temp($node1)\)-$node2\($temp($node2)\)?"

set delta [expr $energy($node1,ENER)-$energy($node2,ENER)]
set delta [expr $delta/(1/($k*$temp($node2)- 1/($k*$temp($node1))))]
set probability [expr exp(-$delta)]
puts "$energy($node1,ENER) $energy($node2,ENER) \nprobability to exchange [expr exp(-$delta)]"

if {$probability >= 1 } { expr 1 } else { expr [prob $probability]  }
}

# production steps 
TERMD begin production

TERMD run 5 TERMDstep exchange temp when {
set k 1.98e-3
puts "exchanging $node1\($temp($node1)\)-$node2\($temp($node2)\)?"

set delta [expr $energy($node1,ENER)-$energy($node2,ENER)]
set delta [expr $delta/(1/($k*$temp($node2)- 1/($k*$temp($node1))))]
set probability [expr exp(-$delta)]
puts "$energy($node1,ENER) $energy($node2,ENER) \nprobability to exchange [expr exp(-$delta)]"

if {$probability >= 1 } { expr 1 } else { expr [prob $probability]  }
}

TERMD destroy
exit
# Alternative way to do replica exchange simulation: change dimension variables instead.
# each node run simulation with constant conformation but exchanging variables
#repeat 100 {
#TERMD do {
#   setSolvent IMM1_memb {temp $temp}
#   redo TERMDstep { temp $temp fname $step numsteps 1000}
#   getEnergy
#} 
#TERMD exchange temp {}
#TERMD exchange com {}
#} 

